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In immiscible two-phase flows, contact line denotes the intersection of the fluid-fluid in- 
terface with the solid wall. When one fluid displaces the other, the contact line moves 
along the wall. A classical problem in continuum hydrodynamics is the incompatibility 
between the moving contact line and the no-slip boundary condition, as the latter leads 
to a non-integrable singularity. The recently discovered generalized Navier boundary con- 
dition (GNBC) offers an alternative to the no-slip boundary condition which can resolve 
the moving contact line conundrum. We present a variational derivation of the GNBC 
through the principle of minimum energy dissipation (entropy production), as formu- 
lated by Onsager for small perturbations away from the equilibrium. Through numerical 
implementation of a continuum hydrodynamic model, it is demonstrated that the GNBC 
can quantitatively reproduce the moving contact line slip velocity profiles obtained from 
molecular dynamics simulations. In particular, the transition from complete slip at the 
moving contact line to near-zero slip far away is shown to be governed by a power-law 
partial slip regime, extending to mesoscopic length scales. The sharp (fluid-fluid) inter- 
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face limit of the hydro dynamic model, together with some general implications of slip 
versus no-slip, are discussed. 



1. Introduction 

The no-slip boundary condition states that there can be no relative motion at the 
fluid-solid interface (Batchelor 1991). It is generally regarded as a cornerstone in con- 
tinuum hydrodynamics, owing to its proven applicability in diverse fluid-flow problems. 
However, decades ago it was discovered that in immiscible two-phase flows, the moving 
contact line (MCL), defined as the intersection of the fluid-fluid interface with the solid 
wall, is incompatible with the no-slip boundary condition (Moffatt 1964; Hua & Scriven 
1971; Dussan & Davis 1974; Dussan 1976; Dussan 1979; de Gennes 1985). As shown 
by Dussan & Davis (1974), under usual hydrodynamic assumptions, viz. incompressible 
Newtonian fluids, no-slip boundary condition, and smooth, rigid, solid walls, there is 
a velocity discontinuity at the MCL, and the tangential force exerted by the fluids on 
the solid bounding surface in the vicinity of the MCL is infinite. This is the well-known 
contact-line singularity. In the past two decades, it was shown through molecular dy- 
namics (MD) simulations that near-complete slip indeed occurs at the MCL (Koplik, 
Banavar & Willemsen 1988; Koplik, Banavar & Willemsen 1989; Thompson & Robbins 
1989; Thompson, Brinckerhoff & Robbins 1993). This finding presented a conundrum 
for classical hydrodynamics, due to a lack of viable alternatives apart from ad hoc fixes. 
Furthermore, in the absence of a viable boundary condition which can reproduce the MD 
results, accurate continuum description of immiscible fiows in the micro- or nanoscales 
remained an elusive goal. 

Through analysis of extensive MD data, it was recently discovered that there is indeed 
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a differential boundary condition, denoted the generalized Navier boundary condition 

(GNBC), which resolves the MCL conundrum (Qian, Wang & Sheng 2003). Here we 

show that the GNBC can be derived variationally from the principle of minimum energy 

dissipation (Onsager 1931a and 1931b), and its implementation through the use of a 

Cahn-Hilliard (CH) free energy functional (Cahn & Hilliard 1958) leads to quantitative 

predictions in excellent agreement with MD simulation results. In what follows, the MCL 

problem is briefly recapitulated in Sec. 2. The variational derivation of the GNBC in Sec. 

3 is followed by a numerical demonstration of its consequences in Sec. 4. It is shown that 

the transition from near-complete slip at the MCL to near-zero slip (no-slip) far away 

from the MCL is not confined to a molecular-sized region around the MCL. Instead, the 

transition follows a power-law profile of partial slipping, extending to mesoscopic scales 

(Qian, Wang & Sheng 2004). The sharp/diffuse (fiuid-fluid) interface limits of our theory 

and their associated issues are described in Sec. 5. In Sec. 6 we discuss some general 

implications of replacing the no-slip boundary condition, which can be regarded as an 

approximation to the GNBC in single-phase flow regions, by the more accurate GNBC. 

In particular, it is argued that slip and partial slip boundary conditions offer the prospect 

of nanoscale interface engineering to "tune" the amount of slipping. 

2. Recapitulation of the moving contact line problem 

Consider an immiscible two-phase flow where one fluid displaces the other (see flgure 
1). If the no-slip boundary condition is applied along the wall, it can be shown that the 
tangential viscous stress varies as r]V/x, where r] is the viscosity, V is the wall speed in the 
reference frame where the fluid-fluid interface is stationary, and x is the distance along the 
wall away from the MCL. This variation leads to diverging stress as a; ^ 0. In particular, 
this stress divergence is non-integrable and implies infinite (viscous) dissipation (Dussan 
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& Davis 1974). Over the years there have been numerous models and proposals aiming 
to resolve this problem. For example, there have been the kinetic adsorption/desorption 
model by Blake & Haynes (1969), the slip models by Hocking (1977), Huh & Mason 
(1977), and Zhou & Sheng (1990), and the diffuse-interface models by Seppecher (1996), 
Jacqmin (2000), Chen, Jasnow & Vinals (2000), Pismen & Pomeau (2000), and Briant 
& Yeomans (2004). 




Figure 1. When one fluid displaces another immiscible fluid, the contact line is moving relative 
to the solid wall. Due to the contact-line motion, the dynamic contact angle dd deviates from 
the static contact angle 9s which is determined by the Young equation 7cosSs -|- 72 =71. 

In another approach to the MCL problem. Cox (1986) carried out an asymptotic 
analysis and found, to the leading order in the capillary number Ca = t]V/j, a relation 
for the dependence of the apparent contact angle (the angle of the fluid-fluid interface 
at some mesoscopic distance away from the fluid-solid interface) on the microscopic 
contact angle, the capillary number, and the distance from the contact line over which 
slip occurs. In this approach, the details of the slip flow in the inner region around 
the MCL are absorbed into a model-dependent constant in the asymptotic relation. 
In particular, it has been shown that the asymptotic behavior is independent of the 
microscopic boundary condition(s) (Dussan 1976). But this conclusion only deepens the 
mystery of what happens at the contact line. 
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In the past two decades, MD simulations have shown that near-complete slip indeed 

occurs at the MCL (Koplik, Banavar & Willemsen 1988; Koplik, Banavar & Willem- 

sen 1989; Thompson & Robbins 1989; Thompson, Brinckerhoff & Robbins 1993). Small 

amount of partial slipping was also observed in single-phase flows (Thompson & Robbins 

1990; Thompson & Troian 1997; Barrat & Bocquet 1999a; Cieplak, Koplik & Banavar 

2001). Such partial slip can be accounted for by the Navier boundary condition (NBC), 

proposed nearly two centuries ago by Navier (1823): t;*'*^' = ^^7, where u**''^ is the slip 

velocity at the surface, measured relative to the (moving) wall, Ig is the slip length, and 7 

is the shear rate at the surface. The small value of Is explains why the NBC is practically 

indistinguishable from the no-slip boundary condition in single-phase macroscopic flows. 

While the NBC can account for the small amount of slip in high shear-rate single- 
phase flows, it fails to account, by an order of magnitude, for the near-complete slip 
at the MCL (Thompson & Robbins 1989; Thompson, Brinckerhoff & Robbins 1993). 
Recently, the MD-continuum hybrid simulation has been applied as a tool to investigate 
this problem by Hadjiconstantinou (1999), and by Ren & E (2005a). But such approaches 
leave unresolved the problem of MCL boundary condition. Lack of a viable boundary 
condition implies that accurate MCL description (and hence immiscible two-phase flows) 
can only be attained through MD simulations, in systems far too small compared with 
most of the experimentally achievable samples. 

The recent discovery of the CNBC (Qian, Wang & Sheng 2003) resolves the MCL 
conundrum. The CNBC states that the slip velocity is proportional to the total tangential 
stress — the sum of the viscous stress and the uncompensated Young stress; the latter 
arises from the deviation of the fluid-fluid interface from its static conflguration. Here 
we show the CNBC to be derivable from the principle of minimum energy dissipation. 
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Its form is hence uniquely necessitated by the thermodynamics of immiscible two-phase 
flows. 

3. Variational derivation of the moving contact line hydrodynamics 

There is a minimum dissipation theorem for incompressible single-phase flows. Accord- 
ing to Batcliclor (1991) (attributable to Hclmlioltz), "the rate of dissipation in the flow 
in a given region with negligible inertia forces is less than that in any other solenoidal 
velocity distribution (of zero divergence) in the same region with the same values of the 
velocity at all points of the boundary of the region." The variational principles involving 
dissipation have been further developed in the works of Rayleigh (1873), Onsager (1931a), 
Onsager (1931b), Edwards & Freed (1974), and Doi (1983). In mathematical terms, if 
we ignore inertia forces for the moment (this can be added in at the end, see equation 
(3.36)) and let the variables ai,- ■ -^an describe the displacement from thermodynamic 
equilibrium, di, • • •, a„ being the corresponding rates and F{a\, ■ ■ ■, an) the free energy, 
then for a set of simultaneous irreversible processes governed by 



where the coefiicients pij are introduced through the linear relations between the rates 
{di} and the "forces" {—dF/dai}, a variational principle (Onsager 1931a and 1931b) 
can be formulated, namely 




(3.1) 




(3.2) 



where ^{d,Q.) is the dissipation- function deflned by 




(3.3) 
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and F{a, a) is the rate of change of the free energy: 

i=l * 

Microscopic reversibihty requires the reciprocal relations ptj = pjt. We note that the 
variation in equation (3.2) should be taken with respect to the rates {di} for prescribed 
{a,}, and the extremum given by equation (3.2) is always a minimum because ^{a, d) is 
quadratic in {d,} and positive definite (Onsager 1931b). The dissipation-function $ as 
given by equation (3.3) is noted to be half the rate of energy dissipation, owing to the 
required consistency between the "force balance" condition (3.1), the minimum condition 
(3.2), and the fact that the dissipative forces are linear in rates. It can be shown that the 
principle of minimum energy dissipation yields the most probable course of an irreversible 
process, provided the displacements from thermodynamic equilibrium are small (Onsager 
1931b, Onsager & Machlup 1953). We note that the variational principle presented here 
is a special form of the principle of minimum energy dissipation formulated by Onsager. 
First, the dissipation- function $ defined here is different from that by Onsager by a 
factor of T, the temperature which is assumed to be uniform in the fluids. Second, the 
rate of change of the free energy F equals —TS+W, where S and W denote the entropy 
and work, respectively (Doi 1983). This point will be elaborated after equation (3.24). 
In Appendix A we use a single-variable case to illustrate the underlying physics of the 
principle of minimum energy dissipation. 

When applied to a single-phase flow conflned by solid surfaces, the variational principle 
in equation (3.2) becomes 

5$(d,d)=^— fci = 0. (3.5) 

Here the rates {di} correspond to the velocity field v(r) and the F{a, d) term deflned by 
equation (3.4) drops out in the variation because there is no such free energy F({aj}) in 
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single-phase flows. As the dissipation- function $(d, d) equals half the rate of energy dis- 
sipation, equation (3.5) leads directly to the minimum dissipation theorem by Helmholtz 
(Batchelor 1991). That is, once the values of the velocity are prescribed at the solid 
surfaces, the rate of viscous dissipation is minimized by the solution of the Stokes equa- 
tion. Physically, when fluid slipping occurs at the solid surface, there is dissipation at 
the fluid-solid interface as well. Below we show that to minimize the total rate of energy 
dissipation Ri [v] (= 2^{a, a.)), the incompressible flow has to satisfy the Stokes equation 
and the NBC simultaneously. Mathematically, this means that the Stokes equation 



-Vp + 7?V^v = 0, 



and the NBC 



pvfP = = -vidnVr + drVn), 

can be derived by minimizing the functional 

'V 



I 



dS 



(3.6) 



(3.7) 



(3.8) 



with respect to the velocity distribution. Here p is the pressure which plays the role of the 
Lagrange multiplier for the incompressibility condition V • v = 0, is the shear viscosity, 
P is the slip coefficient, subscript n denotes the outward surface normal, subscript t 
denotes the direction tangential to the surface, vf^'' is the slip velocity, defined as the 
tangential fluid velocity at the solid surface, measured relative to the (moving) wall, 
al^'' is the nr component of the Newtonian viscous stress tensor, and / dS denotes the 
integration over the solid surface. The functional i?i[v] is 23>(q:, a) for single-phase flows, 
measuring the total rate of dissipation due to viscosity in the bulk and slipping at the 
solid surface. We write the single-phase flow dissipation Ri as the sum of Ry due to 
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viscosity and Rs due to slipping: Ri = Ry + Rs, with 



Rv [v] = J dr [^{diVj + djVifj (3.9) 



and 

R,[v] = j ds[p{vfpf'\ . (3.10) 

While Rs may appear to be very different from R^ in form, yet in reality the two terms are 
very similar if it is reahzed that v^''^^ is just the tangential velocity difference/differential 
at the fluid-solid interface. Hence (3 is just the viscosity divided by a length scale Ig, 
defined as the slip length. Associated with the variation of the velocity field v(r) — > 
v(r) + 5v(r), the change in Ry is given by 

SRy = —2r] J dr [dj {djVi + diVj) dvi] +2ri J dS [dnVrSvr + drVnSvr] , (3-11) 

and that in Rg given by 

SRs = 2p j dS [vfP5vr] . (3.12) 

Imposing the incompressibility condition diVi = by the use of a Lagrange multiplier 
a(r) leads to one more term J dr [adiVi], whose variation is given by 

5 dr [adiVi] = — J dr [diaSvi] . (3.13) 

Here the boundary condition Vn = has been used at the solid surface, where only 
5vt is allowed. From equations (3.11), (3.12), and (3.13), we obtain the Euler-Lagrange 
equations 

-2'qdj {djVi + diVj) - dia = (3.14) 

in the bulk and 

2tj {dnVr + drVn) + 2PvfP = (3.15) 

at the surface. Note that equation (3.14) is identical to the Stokes equation (3.6) with 
djVj = and a = —2p, and equation (3.15) reduces to the NBC (3.7). That equations 
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(3.6) and (3.7) can both be derived variationally from the minimization of -Ri is a gen- 
eralization of the minimum dissipation theorem by taking into account fluid sHpping at 
the solid surface. 

It should be emphasized that while Rs arises from the assumption of fluid-solid in- 
terface slipping, there is no speciflcation of how much sHpping there should be. In other 
words, even an infinitesimal amount of interface slipping would lead to equations (3.6) 
and (3.7). In particular, the no-shp boundary condition is obtained in the limit of /3 ^ oo, 
corresponding to a vanishing sHp length Ig = r]/f3. In the other limit, /? — > 0, we would 
have only the first term on the right-hand side of equation (3.8), and (J^t^ = on the 
boundary. Thus the NBC interpolates between the zero tangential viscous stress limit 
and the no-slip limit. 

To generalize the functional ^{a,d) = ^-Ri[v] from single-phase to immiscible two- 
phase flows, it is recognized that a free energy functional is required to stabilize the 
interface separating the two immiscible fluids. Hence the introduction of a Landau free 
energy functional J^[(?ii(r)] is a necessity, presumably of the form (Cahn & Hilliard 1958; 
Bray 1994) 



where the potential /(</>) has a double-well structure. Here the phase field 0(r) measures 
the (conserved) composition locally defined by ^ = (p2 — Pi)/{p2 + pi), with pi and p2 
being the number densities of the two fluid species. We also introduce the interfacial free 
energy per unit area at the fluid-solid interface, 7/s(0), which is a function of the local 
composition. Two quantities and L can be deflned from the variation of the total free 
energy 




(3.16) 
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that is, 

6 1 J-[</)] + jdS l^fM] } = / I'^^'^] + / [^'^'^] ' (3.17) 
in which ^ = 5T /54> = —KV'^cj) + df{(j))/d(j) by definition is the chemical potential 
and L, given by Kdn^ + d'yfs{^)/d(l), is the corresponding quantity at the solid surface. 
Minimizing the total free energy with respect to 4> yields the equilibrium conditions 
H = C in the bulk and L = at the surface, C being a constant acting as the Lagrange 
multiplier for the conservation of ^. It will be shown later that L = leads to the Young 
equation for the static contact angle. 

From the equilibrium conditions derived above, we see that deviations from the two- 
phase equilibrium may be measured by the "forces" V/x in the bulk and L at the fluid- 
solid interface. For small perturbations away from the equilibrium, the additional rate of 
dissipation arises from system responses that are linear in V/x and L. Such responses 
are described by the diffusive current J in the bulk and the material time derivative of 
(p at the solid surface, i.e., </> = d<f)/dt + Vrdr<i>- The conservation of (j) means that the 
diffusive current and the material time derivative of cj) satisfy the continuity equation 

^=^+wV<i) = -V -3. (3.18) 

But at the fluid-solid interface, diffusive transport normal to the interface is possible 
{dnJn 7^ in general), hence the interfacial (f) is not conserved. In anticipation of the 
relevant dynamics governing the conserved (j) in the bulk and nonconserved interfacial (j), 
J in the bulk and (f) at the solid surface are hence the two additional rates associated 
with the coexistence of two phases. The additional rate of dissipation i?^ due to the dis- 
placement from the two-phase equilibrium may be constructed as a functional quadratic 
in the rates, = Rd + Rr, where 

= rfr ^ , (3.19) 
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1 dS 






r 



(3.20) 

with M and F introduced as two phenomenological parameters. Combining with Ri 
in equation (3.8), we obtain the dissipation-function ^{a,a) for immiscible two-phase 
flows: 



{diVj + djVi 



dS 



f3 



dr 



2M 



dS 



2T 



, (3.21) 



which equals half the total rate of energy dissipation R2 in two-phase flows, i.e., 

$ = --R2 = 2("^l ~'~ ^'t'^ ~ 2^^^ + Rs + Rd + Rr), 



(3.22) 



Now the viscosity 77 and slip coeSicient (3 in equation (3.21), respectively, are understood 
to take on their respective values for the two immiscible fluids on two sides of the interface. 
Note that the right-hand side of equation (3.21) consists of four terms, contributed by the 
four physically distinct sources of dissipation — the shear viscosity in the bulk, the fluid 
slipping at the solid surface, the composition diffusion in the bulk, and the composition 
relaxation at the solid surface. In addition, each term that contributes to $ is positive 
deflnite and quadratic in a rate that arises from the displacement from the equilibrium 
and accounts for a particular source of dissipation. This quadratic dependence follows 
the general rule governing entropy production in a thermodynamic process (Landau & 
Lifshitz 1997), it directly arises from the linear response to a small perturbation away 
from the equilibrium. 

The rate of change of the free energy, i.e., F{a,a) in equation (3.4), may be written 



as 



dv 



d(t> 



dS 



"dt 



(3.23) 



in accordance with the variation of the total free energy in equation (3.17). Substitut- 
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ing d(j)/dt = <j) — V ■ W4> into equation (3.23) and using J dr fxcp = J dr [— /iV • J] = 
/ dr [V/i • J] with J rfr [V • (/uJ)] = J dS [/iJ„] = because of the impermeabihty condi- 
tion J„ = at the sohd surface, we obtain 



F = y dr [V/U • J - MV • V(A] + J dS L{4> - Vrdrcj)) 



(3.24) 



Note that the laws of thermodynamics require F = —TS-\-W, where S and W denote the 
entropy and work, respectively. Here F is the free energy associated with the composition 
field (j), consequently the entropy part —TS must arise from the composition diffusion 
(in the bulk) and relaxation (at the fluid-solid interface) while the work rate W is due 
to the work done by the flow to the fluid-fluid interface. That is. 



-TS = j dr[ViJ.-3]+ j 



dS 



L<j) 



(3.25) 



and 



W = j dr[-v (mV(A)] + j dS [-Vr{Ldr^)] . (3.26) 

It will be seen that /iV0 and Ldrcj) are the "elastic" force/stress exerted by the interface 
to the flow. It is clear that in the steady state F = because the work is fully transformed 
into entropy, i.e., W = TS. This should be obvious since in the steady state the diffusive 
transport of the fluid-fluid interface is balanced by its kinematic transport by the flow, 
hence its free energy is invariant in the course of time. 

Therefore, for immiscible two-phase flows the variational principle in equation (3.2) 
may be expressed by using the functional 



$ + F = / dr {diVj + djVif\ +fdS 



Jdr 



2M 



JdS 



2T 



+ 



Jdr[\/ii-3 - fiv +JdS L{4> - Vrdr4>) 



(3.27) 

Based on equation (3.27), a hydrodynamic model for the contact-line motion can be 
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derived by minimizing $ + f with respect to the rates {v, J, (j)}, supplemented with the 
incompressibihty condition V • v = 0. 

As $ is quadratic in (J, cj)) and F is Unear in (J, </)), the Euler-Lagrange equation with 
respect to J is given by 

J = -MV/u, (3.28) 

where the parameter M introduced in equation (3.19) is seen to have the meaning of a 
mobiUty coefficient. Substituting equation (3.28) into the continuity equation (3.18) for 
(f) yields the anticipated advection-diffusion equation 

^=^+vV(A = -V-J = MVV- (3.29) 

Similarly, the corresponding Euler-Lagrange equation for minimizing $ + F with respect 
to (j) at the solid surface is 

= ^ + Vrdr^ = -TL{4>). (3.30) 

That is, at the fluid-solid interface, the relaxation dynamics of the interfacial </> is linear 
in L{(f)), i.e., AUen-Cahn dynamics for nonconserved quantities (Bray 1994). 
Now we show that the Stokes equation with the capillary force, 

-Vp + r]V\ + ijlV(I) = Q, (3.31) 

and the GNBC with the uncompensated Young stress, 

/3(<^)t;f f = -ll{dnVr + drVn) + L{(l>)dr<l>, (3.32) 

can be obtained by minimizing $ + F with respect to the fluid velocity. Here ;uV0 is 
the capillary force density (Chella & Vinals 1996; Qian, Wang & Sheng 2003), (3{<p) is 
the slip coefficient which may locally depend on the composition, and L{(j))dr4> is the 
uncompensated Young stress which vanishes in equilibrium (Qian, Wang & Sheng 2003). 
From equations (3.21), (3.22), (3.24) and (3.27), we see that the dependence of on 
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the velocity comes from = ^{Ry+Rs) in $ and / dr [— /iv • W4>\ + J dS [—Lvrdrcf)] in 

F. Consider a variation of the velocity field v(r) — > v(r) + 5v(r). The associated changes 

in Ry and Rg are already given by equations (3.11) and (3.12), and those in F are given 

by 

6F = j dr [-/iV^ ■5v]+ j dS [~Ldr(l>Svr] = j dv [- ijLdi(t)5vi] + j dS [-Ldr(j>Svr] • 

(3.33) 

Combining equations (3.11), (3.12), (3.13), and (3.33), we obtain the Euler-Lagrange 
equations 

-r]dj {djVi + diVj) - ^dia - ^dicf) = (3.34) 

in the bulk and 

r] [dnVr + drVn) + 0vfP - L{(j))dr<f> = (3.35) 

at the surface. Note that equation (3.34) is identical to the Stokes equation (3.31) with 
djVj = and a = —2p, and equation (3.35) reduces to the GNBC (3.32). An important 
point of this derivation is that the uncompensated Young stress at the boundary (last 
term on the left side of equation (3.35)) must accompany the capillary force density in 
the bulk (last term on the left side of equation (3.34)), both being the "elastic" interfacial 
force. Hence the uncompensated Young stress is simply the manifestation of the fluid-fluid 
interfacial tension at the solid boundary. 

Once the free energies J^[(j>] and 7/s(?5') are fixed, the contact-line motion (in the regime 
of small Reynolds number) is fully determined by equations (3.29), (3.30), (3.31), and 
(3.32), supplemented by the incompressibility condition V-v = and the impermeability 
conditions Vn = and 9„/U = at the solid surface (Qian, Wang & Sheng 2003). The 
Stokes equation can be readily generalized to the Navier-Stokes equation 



, „ ' 
^ + (v.V)v 



= -Vp + ryV^v + yN4>, (3.36) 
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by including the inertia forces, where p is the mass density. Together, the Navier-Stokes 
equation (3.36), the GNBC (3.32), the advection-diffusion equation (3.29), and equation 
(3.30) for the relaxation of interfacial (f), form a consistent hydrodynamic model for the 
contact-line motion in immiscible two-phase flows, first presented by Qian, Wang & Sheng 
(2003). It is now clear that our model is necessitated by more general considerations. 



4. Comparison between MD and continuum results 

To demonstrate the physical validity of our model, numerical solutions have been ob- 
tained for direct comparison to the MD velocity and interfacial profiles. For this purpose, 
we make use of the CH free energy functional (Cahn & Hilliard 1958) 



-^Cff[</>(r) 



dr 



(4.1) 



to fix the form of f{^) for J^[^] in equation (3.16). Here K, r, and u are material parame- 



ters that can be determined from the interfacial thickness ^ — y'K/r, the interfacial ten- 
sion 7 = 2\/2r^^/3u, and the two homogeneous equilibrium phases ^± = ^^prju = ±1. 

The two coupled equations of motion are the advection-diffusion equation for the phase 
field </)(r) and the Navier-Stokes equation in the presence of the capillary force density: 



^+vV</. = MVV 



(4.2) 



9v 

at 



-h (v • V) V 



-Vp + V-CT^+^V^ + fe, (4.3) 

together with the incompressibility condition V-v = 0. Here M is the mobility coefficient, 
[i = 5TchI^4> is the chemical potential derived from the CH free energy functional Tch^ 
p is the mass density of the fluid, p is the pressure, cr" = [(Vv) -|- (Vv)-^] is the 
Newtonian viscous stress tensor with 77 being the viscosity, jiVcj) is the capillary force 
density, and fe is the external force. The boundary conditions at the solid surface are the 
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impermeability conditions ^n^l = 0,Vn = 0, the relaxational equation for surface (j): 



Here r denotes the direction tangent to the solid surface, n denotes the outward sur- 
face normal, F is a positive phenomenological parameter, L{(j)) = Kdncj) + d'yfs{(t>)/d(j) 
with 7/s(?!') being the fluid-solid interfacial free energy per unit area, /?(</)) is the slip 
coefficient which may locally depend on the local surface composition 0, and L((p)dr<j> 
is the uncompensated Young stress. We use 7/s(?!>) = (A7/s/2) sin(7r(/)/2) which is a 
smooth interpolation from 7/s(0_) = — A7/s/2 to 7/s(<?i+) = A7/s/2. According to the 
Young equation for the static contact angle Os'- ^fs{<p+) + jcosOg = 7/s((/)_), we have 
A7/S = -7C0S^s. 

To interpret physically the second term on the right-hand side of equation (4.5), let us 

consider a fluid-fluid interface that intersects the planar solid surface z = with a contact 
angle 9 relative to the x axis. For simplicity we assume the contact line to be straight 
and parallel to the y axis, and hence r = a;. If the fluid-fluid interface is gently curved, 
then J^^_^dx[{Kdn4')dx4'] — Ji„i d4>{K 8^4') (^os 9 , where J^^^dx denotes the integration 
across the fluid-fluid interface along x and dm means taking spatial derivative along the 
fluid-fluid interface normal m, with dn^ ~ dm^cosO. As J^^^d(j){Kdm4') = 7) we have 
J.^j dx [{Kdn(l>)dx4>\ = 7cos^. In the equilibrium, L{(j)) = 0, and hence 




FL(</>), 



(4.4) 



and the GNBC in continuum differential form: 



Pi^K^ = -vidnVr + drVn) + L{cj))dr(l>. 



(4.5) 




vanishes. This leads directly to the Young equation 



7 cos 6*^-^7/^ ((/>+) -7/^ ((/)_) =0, 



(4.6) 
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where 6s is the static contact angle, 7 cos 9s comes from J.^^ dx [(i^9„^)9x^] , and 7/s (</>+)— 
7/s(^_) from J.^j dx [{d'jfs/d(p)dx4']- When the fluids are in motion, integrating the un- 
compensated Young stress L{(l))dx(j> across the fluid- fluid interface along x yields 

/ dx[L{(t))da:(f)]=jcosed + ^fsi<P+)-7fsi(l>-)=l{cosed-coses), (4.7) 

J int 

where 64 is the dynamic contact angle. Equation (4.7) implies that the uncompensated 
Young stress arises from the deviation of the fluid-fluid interface from its static config- 
uration. However, it must be pointed out that here the contact angle is the so-called 
"microscopic contact angle." It difi^ers from the "apparent contact angle" when the con- 
tact line is in motion. In particular, the apparent contact angle can change rather sharply 
when the contact line begins to move. But 9d can only vary linearly and relatively slowly 
with velocity. 

Prom equation (4.7) and the fact that for moderate flow rates the field (j) for gently 
curved (fiuid-fluid) interface relaxes to essentially the local stationary structure, we can 
write 

L{(t>)dx(t> = jcos0df{x) + dx^fs, (4.8) 
where f{x) is a function which peaks at xcl, the contact line center, and / dxf{x) = 1. 
In particular, when the static contact angle is 90° and 7/s is a constant, f{x) can be well 
approximated by (3/4-\/2^) cosh~^[(a; — xcl)/V^^], derived from the stationary solution 
= (l>+ tanh[(a;— a;ci)/-\/2^]. Thus the uncompensated Young stress is a peaked function 
centered at the contact line (note that 9x7/s is also peaked at the interfacial region). 
While the above expressions give a physical interpretation to the uncompensated Young 
stress, it is noted that the presence of cos^d makes the expression awkward for the 
purpose of calculation (as 6d should be the outcome of the calculation). Thus in actual 
computations the formulation in terms of (/> is prefered. 
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We have carried out the MD-continuum comparison in such a way that virtually no 

adjustable parameter is involved in the model calculation (Qian, Wang & Sheng 2003). 
This is achieved as follows. There are a total of nine material parameters in our model: 
p, T], (3, ^, 7, \(j>±\, 9s M, and F. Among these parameters, the first seven are directly 
measurable in MD simulations. As for M and F, their values are fixed through an op- 
timized MD-continuum comparison: one flow field from MD simulation is best fitted by 
that from hydrodynamic model calculation with optimized M and F values, although in 
our case the fitting is not very sensitive to the values of M and F. That is, fitting can 
be almost equally good if M and F deviate from the optimal values. It will be shown in 
Sec. 5 that such insensitivity is not an accident, i.e., as long as the values of M and F 
are in the right range for the hydrodynamic model to reach the sharp interface limit, the 
continuum predictions should not be sensitive to these parameter values. 

Once all the parameter values are determined, predictions from our hydrodynamic 
model can be readily compared to the results from a series of MD simulations with dif- 
ferent external conditions. The overall agreement is excellent, thus demonstrating the va- 
lidity of the GNBC and the hydrodynamic model. We emphasize that the MD-continuum 
agreement has been achieved from the molecular-scale vicinity of the contact line (Qian, 
Wang & Sheng 2003) to far fields at the large scale (Qian, Wang & Sheng 2004). 

MD simulations have been carried out for immiscible two-phase flows in Couette ge- 
ometry (see figure 2) (Qian, Wang & Sheng 2003). Two immiscible fiuids were confined 
between two planar solid walls parallel to the xy plane, with the fluid-solid interfaces 
defined at 2 = and H. The Couette flow was generated by moving the top and bottom 
walls at a constant speed V in the ±x directions, respectively. Periodic boundary condi- 
tions were imposed along the x and y directions. Technical details of our MD simulations 
may be found from Qian, Wang & Sheng (2003) and Qian, Wang & Sheng (2004). Two 
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cases were considered. In the symmetric case, the static contact angle 9s is 90° and the 
fluid-fluid interface is flat, paraflel to the yz plane. In the asymmetric case, the static con- 
tact angle 9s is 64°, and the fluid- fluid interface is curved in the xz plane. Steady-state 
velocity and interfacial profiles were obtained from time average over lO^r or longer, 



where r is the atomic time scale ^Jma'^/e, with e and a being the energy and length 
scales in the Lennard-Jones potential for fluid molecules, and m the fluid molecular mass. 
Throughout the remainder of this paper, all physical quantities are given in terms of the 
Lennard-Jones reduced units (defined in terms of e, a, and m). 

(a) , 

H 



fluid 1 



fluid 



flluidl 




synunetric case 
(b) 



fluid 

asymmetric case 




Figure 2. Schematic of simulation geometry, (a) Static configurations in the symmetric and 
asymmetric cases. Fluid 2 is sandwiched by fluid 1 due to the periodic boundary condition along 
the X direction, (b) Dynamic configuration in the symmetric case. 

Figure 3 shows the MD and continuum velocity fields for a symmetric case of Couette 
flow, figure 4 shows those fields for an asymmetric case of Couette flow, and figure 5 
shows the MD and continuum fluid-fluid interface profiles for the above two cases, from 
which the difference between the "apparent" advancing and receding contact angles is 
clearly seen. 

In order to further veriiy that the hydrodynamic model is local and the parameter 
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Figure 3. Comparison of the MD (symbols) and continuum (lines) velocity profiles {vx{x) and 
Vz{x) at different z levels) for a symmetric case of immiscible Couette flow {V = 0.25(e/m)^''^ 
and H = 13.6cr). The profiles are symmetric about the center plane z = H/2, hence only the 
lower half is shown at z = 0.425cr (circles and solid lines), 2.125cr (squares and dashed lines), 
3.825cr (diamonds and dotted line), and 5.525cr (triangles and dot-dashed lines). 

values are local properties, hence applicable to different flow geometries, we have carried 
out MD and continuum simulations for immiscible Poiseuille flows (Qian, Wang & Sheng 
2003). We find that the hydrodynamic model with the same set of parameters is capable 
of reproducing the MD velocity and interfacial profiles, shown in figure 6. Similar to 
what's observed in Couette flows, here the slip is near-complete at the MCL, i.e., ^ 
and |t'j'*''| « V (the wall speed), while far away from the contact line, the flow field is not 
perturbed by the fluid-fluid interface and the single-fluid unidirectional Poiseuille flow is 
recovered. 

The parameter values used for the MD-continuum comparison have been given by 
Qian, Wang & Sheng (2003). We emphasize that the overall agreement is excellent in all 
cases, therefore the validity of the GNBC and the hydrodynamic model is well affirmed. 

It should be mentioned that our MD results also show some fluid-fluid interfacial 
structure which is outside the realm of our continuum model. In particular, the MD data 




Figure 4. Comparison of the MD (symbols) and continuum (lines) velocity profiles {vx{x) and 
Vz{x) at different z levels) for an asymmetric case of immiscible Couette flow {V = 0.2(e/m)^^^ 
and H = 13.6cr), shown at z = 0.425cr (circles and solid lines), 2.975cr (squares and long-dashed 
lines), 5.525cr (diamonds and dotted line), 8.075cr (up-triangles and dot-dashed lines), 10.625cr 
(down-triangles and dashed lines), 13.175cr (left-triangles and solid lines). Although the solid 
linos arc used to denote two different z levels, for each solid line, whether it should be compared 
to circles or left-triangles is self-evident. 
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Figure 5. Comparison of the MD (symbols) and continuum (lines) fluid-fluid interface profiles, 
defined by pi = P2 (0 = 0). The circles and dotted line denote the symmetric immiscible Couette 
fiow with V = 0.25(e/m)^''^ and H = 13.6a; the squares and dashed line denote the asymmetric 
immiscible Couette flow with V = 0.2(e/m)^^^ and H = 13.6cr. 
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x/a 

Figure 6. Comparison of the MD (symbols) and continuum (lines) results for an asymmetric 
case of immiscible Poiseuille flow. An external force mgext = 0.05e/<T is applied on each fluid 
molecule in the x direction, and the two walls, separated hy H = 13.6cr, move at a constant 
speed V = 0.51(e/m)^''^ in the —x direction to maintain a stationary steady-state interface. 
Upper panel: Fluid-fluid interface profiles, deflned by pi = p-i (</> = 0). Lower panel: Vxix) at 
different z levels. The profiles are symmetric about the center plane z — H/2, hence only the 
lower half is shown at z = 0.425(T (circles and solid line), 2.125(7 (squares and dashed line), 
3.825a (diamonds and dotted line), and 5.525cr (triangles and dot-dashed line). 

show a fairly significant density drop right at the interface (</> = 0). Physically this is due 
to the repulsive molecular interaction between the two immiscible fluids. This density 
drop has implication on the slip coefficient (3 in that for the interfacial region (3 is no 
longer a simple composition of /3i and (32 on the two sides of the contact line. At the 
same time, MD data also revealed and verified the form and nature of the uncompensated 
Young stress, just as predicted by equations (4.7) and (4.8). Details can be found in Qian, 
Wang & Sheng (2003). 

Another important comparison between MD and continuum hydrodynamics is the 
behavior that interpolates between the MCL, where there is near-complete slip, and far 
away from the contact line, where there is at most a small amount of partial slip. 
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MD simulations have been carried out for immiscible Couette flows in increasingly 

wider channels (Qian, Wang & Sheng 2004). The inset to figure 7 shows the tangential 

velocity profiles at the wall. Immediately next to the MCL, there is a small core region, 

where the slip profiles show a sharp decay within a few Ig (the slip length in single-phase 

flows). As the channel width H increases, a much more gentle variation of the slip profiles 

becomes apparent. In order to reveal the nature of this slow variation, we plot in figure 

7 the same data in the log-log scale. The dashed line has a slope —1, indicating the 1/x 

behavior of the slip profile, where x is the distance from the MCL. Because of the finite H, 

there is always a plateau in each of the single-phase flow regions, where the constant small 

amount of slip is given by Uq'*^ = 21//s/(iJ-|-2/s), which acts as an outer cutoff on the 1/x 

profile. This expression for Vq^^ is simply derived from the NBC and the Navier-Stokes 

equation for uniform shear flow. Our largest MD simulation shows that the i)*''^ oc 1/x 

behavior extends to ~ 50(7 (or ~ 25/s). Obviously, as — > oo and Vq^^ approaches 

(no-slip), the power-law region can be very wide indeed. This large 1/x partial slip region 

indicates that the outer cutoff length scale (e.g., the system dimension) would determine 

the integrated effects, such as the total steady-state dissipation. Actually, in the past the 

similarity solutions of Stokes equation have shown the 1/x stress variation away from 

the MCL (Moffatt 1964; Hua & Scriven 1971). However, to our knowledge the fact that 

the partial slip also exhibits the same spatial dependence was first determined by Qian, 

Wang & Sheng (2004), even though the validity of the NBC has been verified at high 

shear stress (Thompson & Troian 1997; Barrat & Bocquet 1999a; Cieplak, Koplik & 

Banavar 2001). 

The continuum results shown in figure 7 were obtained on a uniform mesh, using the 
same set of parameter values corresponding to the same local properties in all the five 
MD simulations. We have extended the MCL simulations, through continuum hydrody- 
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Figure 7. Log- log plot for the slip profiles. Here v^/V + 1 is the scaled slip velocity at the lower 
fluid-solid interface z = 0, and x/a measures the distance from the MCL in units of a. The wall 
is moving at —V, hence v^/V — means complete slip and Vx/V = —1 means no slip. The 
Vx profiles were obtained for five symmetric cases of Couette fiow, which used different values 
for H but the same value for V {— 0.05-y/ e/m) and also the same parameters for densities and 
interactions. The symbols represent the MD results and the solid lines represent the continuum 
results, obtained for H = 6.8a (black circles and line), H — 13.6cr (red squares and line), 
H = 27.2(T (green diamonds and line), H = 54.4(t (blue up-triangles and line), H — 68a (orange 
down-triangles and line). There are two solid curves for each color, one for the slip profile left 
to the MCL and the other for that right to the MCL. The dashed line has the slope of — 1, 
indicating that the 1/x behavior is approached for increasingly larger H. For H — 68a, the 1/x 
behavior extends from |a;| ~ 12a ~ 6I3 to 50a ~ 25ls, where Is was measured to be 2a. Inset: 
The scaled tangential velocity Vx/V at 2 — 0, plotted as a function of x/a. 

namics, to lower flow rates and much larger systems. For this purpose, we have employed 
the adaptive method based on iterative grid redistribution (Ren & Wang 2000). The 
computational mesh is redistributed following the behavior of the continuum solution so 
that fine molecular resolution is achieved in the interfacial region, while elsewhere a much 
coarser mesh is used to save computational cost. A semi-implicit time stepping scheme 
is also used to speed up the approach to steady state. 
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Figure 8. Log- log plot of the slip profiles. Here Vx/V -I- 1 is the scaled slip velocity at the 
lower fluid-solid interface z = Q, and x/ls measures the distance from the MCL in unit of 
la- The wall is moving at —V, hence Vx/V — means complete slip and Vx/V — —1 means 
no slip. There are three symmetric cases of Couette flow. The black solid line denotes the 
case of H — 326a, V = 0.005-y/ e/m and h = 1.24(j, the red dashed line denotes the case of 
H = 326a, V = 0.0025-^ e/m and h = 1.24a, and the blue dotted line denotes the case of 
H = 326a, V = 0.0025-y/ e/m and Is = 0.62a. The green dot-dashed line has the slope of —1, 
indicating a 1/a; region much wider than that in figure 7. Inset: Universal slip profile. The scaled 
tangential velocity Vx/V at z = is plotted as a function of the scaled coordinate x/ls for the 
three cases. The slip profiles exhibit a partial-slip region as large as hundreds of Is - The relation 
Vx/V = 1/(1 + |x|/2.14Zs) — 1 is also plotted by the orange dot-dashed line, showing an extremely 
good fit. 

The continuum results for three large systems of low flow rates are shown in figure 8. 
The capillary force is verified to be important in the interfacial region only, whereas the 
pressure gradients and viscous forces show a much slower variation. They are balanced 
outside the interfacial region, and hence the flow is governed by the Stokes equation. 
This is expected, because the Reynolds number mpYH/rj w 0.6 for p sa O.S/cr^, V = 



0.005v/eM: H w 300(7, and w 2.0^/^/(7'^. Fi gure 8 shows the slip profiles plotted 
on the log- log scale, which clearly exhibit the 1/x behavior extending from |x| « Qlg 
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to \x\ K, 2701s. The inset to figure 8 shows the scaled tangential velocity profiles at the 

solid surface, from which the existence of a universal slip profile is evident. Physically, 

when H ^ Is, the Stokes flow is governed by only one velocity scale V and one length 

scale Is- Thus universality becomes evident from v^/V plotted as a function oi x/ls- A 

heuristic account of the universal slip profile is as follows. Away from the MCL, the 

viscous shear stress is given by -'arjVx{x)/\x\, where a is a constant ~ 1 and Vx{x) is 

the local tangential velocity. The NBC implies v^'^^ix) = —alsVx{x)/\x\. Combining this 

equation with vf^'P = v^ + V yields v^J^'p{x)/V = 1/(1 + \x\/als). This relation, with 

a K, 2.14 for best fit, agrees with the continuum slip profiles extremely well, as seen in 

the inset to figure 8. 

5. Scaling analysis: sharp or diffusive interface limit 

In this section we look at the MCL problem from two different perspectives. These 
naturally arise from the physical reality that there are two distinct regions in the MCL 
problem: the interfacial region and the rest. In the perspective of the sharp interface limit, 
the problem is looked at from outside the interfacial region. The effect of the interfacial 
region is taken into account only in an integrated sense. In the diffuse interface limit, on 
the other hand, we examine the problem when the interfacial region is physically large. It 
turns out that the two limits can be mathematically analyzed by varying the magnitude 
of M and F. In addition, we also give an abbreviated account of the recent work by Ren 
& E (2005b) on the sharp interface limit, which not only simplifies the mathematics, but 
also expresses the GNBC without the (j) variable, and hence physically more transparent. 

In appendix B we have derived the total rate of energy dissipation in steady state (see 
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equations (B3), (B4), and (B5)), 

R2[^] = Rvbr]+Rs[v]+Rd[v]+Rr[v] 



+ JdS 



^ {Vrdr(t)f 



(5.1) 

This function(al) is useful to the variational analysis of the sharp/diffuse interface limits. 
It is obtained by eliminating the two rates (J, </>) by expressing them in v, for the steady 
state. The Stokes equation (3.31) and the GNBC (3.32) can be derived by minimizing 
-R2[v] with respect to v, supplemented with equations (3.29) and (3.30) for the linear 
dissipative dynamics of </> (see equations (B8) and (B9)). 

The CH diffuse-interface modeling allows diffusive transport through the fluid-fluid 
interface (Seppecher 1996; Jacqmin 2000; Chen, Jasnow & Vinals 2000; Pismen & Pomeau 
2000; Briant & Yeomans 2004). By taking the limit of M ^ and T ^ 0, we obtain the 
sharp interface limit in which advection dominates such that ^ — > and d(p/dt — > — v- V^. 
According to equation (5.1), if 1/M and l/F both approach positive infinity, minimizing 
R2 with respect to v then requires v • V0 ^ in the interfacial region (see the last two 
terms on the right-hand side of equation (5.1)). That means in steady state the flow is 
parallel to the interface. 

The limiting magnitude of v- V0 can be obtained through a scaling analysis. We assume 
that the interfacial thickness ^ is much smaller than the smallest interfacial curvature 
radius, a limit realized for moderate shear rates. Integrating the Stokes equation (3.31) 
across the fluid-fluid interface yields ~ ?7||Vv||, from which wc have ~ riV/l(j)+, 

where V and I are the characteristic velocity and length scales of the flow. Integrating 
equation (3.29) (with d(j)/dt = in steady state) across the fluid-fluid interface yields 
~ -^||V/x||, from which we have \vm\ ~ -M|a*|/^0+) where Vm is the velocity com- 
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ponent along the interfacial normal m. Together, r]V/l(j)+ and \vm\ ~ M|/x|/Z(/)+ 

lead to It^rol/^ ^ Mr]/P(ly^ for the dimensionless interfacial normal velocity. This rela- 
tion indicates that in order to realize the sharp interface limit, ^/Mrj/<j>^ is the length 
scale which must be made small enough compared to I. (This length scale arises from the 
coupling of equations (3.31) and (3.29) (Bray 1994; Jacqmin 2000).) As for the charac- 
teristic length / of the flow field, it is given by the slip length Ig in the inner region close 
to the contact line, where Is is actually the only length scale governing the Stokes flow. 
Far away from the contact line, I becomes the dimension of the conflned system, e.g., 
I ^ H. Typically Ig is much smaller than H, and consequently \vm\ reaches the maximum 
in the immediate vicinity of the MCL. This leads to the first conclusion that the sharp 
interface limit is realized in the bulk when 

This relation ensures the flow to be parallel to the interface in the inner region of fast 
velocity variation and large interfacial curvature. 

In the sharp interface limit, the tangential viscous stress is vanishingly small at the 
MCL because of the vanishing interfacial normal velocity there. Here we consider the 
case of 9s = 90° for simplicity and generalization to 9s 90° is straightforward. It 
follows that in steady state the near-complete slip at the MCL is mainly sustained by 
the uncompensated Yoiuig stress: /3V ^ -y cos 9d/£,, where V is for the maximum slip 
velocity (the wall speed relative to the stationary fluid-fluid interface) and ycos9d/i is 
obtained by considering the integrated uncompensated Young stress in equation (4.7) 
as distributed within ^. The difference between the dynamic and static contact angles, 
= 9s — 9d, is then obtained as A9 ~ ^PV/j = ^Ca/ls, where Is = r}/(3 is the slip length 
and Ca = r]V/'y is the capillary number. Integrating equation (3.30) (with d(p/dt = 0) 
across the fluid-fluid interface along the solid surface yields the tangential velocity across 
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the interface |^^r| Fi^cos^d « TKA6. (For 6s = 90°, the tangential direction r is the 
same as the interfacial normal m, and hence Vr becomes Vm-) Substituting ^Ca/ls 
into \vm\ ~ TKA0, we obtain \vm\ ~ TK^Ca/lg, and hence luml/V' ~ ^K^rj/jls- Using 
7 = 2y/2K<j)'^/3^, we obtain [utoI/'^^ ~ T^^r]/(j)^ls. This leads to the second conclusion 
that the sharp interface limit is realized at the contact line when 

which ensures the interface not to be penetrated by the flow. Note that the two conditions 
(5.2) and (5.3) are independent of the magnitude of V. This is consistent with the linearity 
of the hydrodynamic model (for sufficiently low flow rates) and has been well verified 
numerically. 

In Sec. 4, it has been mentioned that M and F are treated as fitting parameters 
to optimize the MD-continuum comparison. Since the fluid-fluid interface in our MD 
simulations is impenetrable, the hydrodynamic model has to be within the sharp interface 
limit in order to reproduce the MD results. It follows that the values of M and F must 
satisfy the conditions (5.2) and (5.3). Moreover, as long as the sharp interface limit is 
reached, the continuum predictions are not sensitive to the values of M and F. Such 
has indeed been our experience. That means M and F should not be regarded as fitting 
parameters in our MD-continuum comparison because they are simply used to realize the 
sharp interface limit, in accordance with the interface impenetrability conditions (viewed 
outside the interfacial region). 

We have shown that in order to sustain the near-complete slip at the MCL, the dynamic 
contact angle 9d has to deviate from the static angle 9s by ~ ^Ca/ls- A similar 
scaling relation has recently been obtained by Ren & E (2005b) in a general discussion 
on the sharp interface limit of the GNBC. Here we derive a complete expression for 
within the CH phase-field formulation. Consider an interface deformed by the shearing 
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movement of confining walls (see figure 2b). In steady state the tangential viscous stress is 

negligibly small at the MCL where the slip is near-complete (v*''*' w F at the lower fluid- 
solid interface). Therefore the GNBC may be expressed as 0V = —Kdz(j>dx(j)- Suppose 
the deviation of 9^ from 9^ = 90° is very small (this is indeed the case, as we will 
show later). Then ~ —dx(j>cos9d and w K {dx4>f cos, 9 d- Using the interfacial 
profile ^(a;) = ianh {x / \/2^) at 2 = (for a gently deformed interface), we obtain 
13V « K(l)\ cos6'd/2^2_ follows that 

This expression has been quantitatively verifled. Substituting the MD values Ig = 3.8^, 



7 = 5.5e/(T^, r] = IM^/em/a^, and V = ^.I^^JTJm into equation (5.4) yields = 0.044 
(or 2.5°), in excellent agreement with cos^d = 0.0437 obtained from the ^ = locus in 
the continuum solution. Actually there is another measure of cos 9d using the integrated 
Young stress J.^^dxKdn4'dx4> = 7C0s6'd, which produces a value only slightly different 
from that determined by the ^ = locus. We note that in our nanoscale MD simulations. 
Is is typically larger than ^ and 7 > 10r]V, and hence A9 < 0.1. Yet this small angle of 
deviation is necessary for the uncompensated Young stress to sustain the near-complete 
slip of the MCL (Qian, Wang & Sheng 2003). 

In a recent study of the sharp interface limit of the GNBC by Ren & E (2005b), it has 
been shown that the deviation of the dynamic contact angle 9d from the static contact 
angle 9s is proportional to the dimensionless parameter (3*V5/^, which measures the 
relative strength of the frictional force between the fluid and the solid and the interfacial 
force between the fluids. Here /3* is the (average) slip coefflcient in the contact line 
region, depending on /3i and (32 in the two single-phase flow regions as well as the fluid- 
fluid interfacial structure, and 5 is the fluid-fluid interfacial thickness. Numerical results 
have been obtained for the relation between the microscopic contact angle 9d and the 
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apparent contact angle, demonstrating good agreement with the analytical results based 
on matched asymptotic expansions (Cox 1986). 

Prom the above, it follows that in the sharp interface limit the contact angle can be set 
at the value of the static contact angle (for ^ ^ and/or Ca ^ 0), and since outside the 
fluid-fluid interfacial region the NBC is valid, numerically the flow fleld can be calculated 
separately on the two sides of the interface, linked together via the interface transition 
relations (zero normal component of velocity, continuity of tangential component of ve- 
locity, continuity of tangential stress, normal stress difference across the interface being 
balanced by the tensile force proportional to the interface curvature) (Zhou & Sheng 
1990; Ren & E 2005b). However, in such numerical solutions the tangential (viscous) 
stress is necessarily discontinuous at the contact line. This can be clearly seen by consid- 
ering the tangential viscous stresses approached along the three interfaces (two fluid-solid 
and one fluid-fluid) terminating at the contact line. For simplicity, let us consider a fluid- 
fluid interface vertically intersecting the solid surface. Velocity continuity dictates that 
the tangential viscous stress at the solid surface approaches (3iV (or P2V) at the contact 
line in the single-phase flow region left (or right) to the fluid-fluid interface, with (3i 
and /32 being the slip coefflcients in the left and right regions, respectively. However, the 
impenetrability condition at the fluid-fluid interface dictates that along this interface, 
the velocity component parallel to the solid surface vanishes, leading to a vanishing tan- 
gential viscous stress at the solid surface when approached along the fluid-fluid interface 
down to the contact line. Hence there exist three distinct values for the tangential viscous 
stress. The uncompensated Young stress thus enters as the required subsidiary condition 
to complete the picture and make the solutions physically meaningful. In particular, the 
uncompensated Young stress interpolates between the two values of piV and (32V, with 
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a mean value given by {(3i + (32)V/2, if the two fluids are assumed to interact with the 

solid independently (Qian, Wang & Sheng 2003). 

Opposite to the sharp interface limit is the limit of M — > oo and F — > oo. Obviously, in 
this limit minimizing i?2 in equation (5.1) is equivalent to minimizing i?i in equation (3.8) 
because Rd and Rr both vanish regardless of the velocity distribution (see equations (B 3) 
and (B4)). As a consequence, the flow field and the fluid- fluid interface are decoupled: 
the velocity is distributed as if there is only one single phase while the interfacial profile 
approaches the equilibrium one. As the fluid-fluid interface becomes very transparent 
(through difi^usive transport), the contact line loses its usual implications. Indeed, it has 
been shown that as this limit is approached, the stress singularity can be lifted even if 
the no-slip boundary condition is applied (Seppecher 1996; Jacqmin 2000; Chen, Jasnow 
& Vinals 2000; Briant & Yeomans 2004). This result can be made intuitively plausible 
from our variational formulation as follows. 

According to equation (5.1), if the functional i?2[v] is to be minimized subject to /3, 
1/M, and l/F all approaching positive infinity, then the problem is reduced to solving 
the Stokes equation subject to the no-slip boundary condition vf^P = and the interface 
impenetrability condition Vm = 0. This leads to the well-known non-integrable singular- 
ity in viscous dissipation. Therefore, mathematically the contact-line singularity may be 
viewed as resulting from minimizing i?2[v] with /3, 1/M, and l/F — » oo. By removing 
either the /? — > oo constraint (i.e., allowing slipping), or the 1/M, l/F ^ oo constraint 
(i.e., allowing diffusive relaxation), the total dissipation can only decrease from infinity, 
thus regularizing the solution. This is especially the case since the divergence is loga- 
rithmic in nature, i.e., the divergence is only marginal. We should note, however, that 
physically realistic cases correspond to (3, 1/M, and l/F all remain finite, as evidenced 
by MD results. In fact, application of the above considerations to the problem of corner- 
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flow singularity (involving a flow in a corner with one rigid plane sliding over another) 
(Batchelor 1991; Moffatt 1964; Kophk & Banavar 1995) would be equally vahd (Qian & 
Wang 2005). 

It is interesting to note that in either of the two limits discussed above, the rate of 
interfacial dissipation, Rd+Rr, tends to vanish. In the sharp interface limit of M — > and 
r ^ 0, the limiting behaviors of the interfacial normal velocity expressed in equations 
(5.2) and (5.3) make i?d ~ M and Rr ^ F. (Equation (B 3) indicates Rd oc \vm\'^/M while 
equation (5.2) indicates \vm\ oc M, and hence ii^ oc M —> 0. Similarly, Rj. ocT ^ 0.) In 
the opposite limit of M — > oo and F ^ oo, as the interface is penetrated by the flow, 
Rd and Rr simply vanish as 1/M and 1/F, respectively. That the positive deflnite rate 
of interfacial dissipation, Rd + Rr, approaches zero in the two opposite limits implies 
that a maximum should be reached somewhere in between. However, the total rate of 
dissipation i?2 should increase monotonically from the limit of M oo and F ^ oo 
to that of M — > and F — > 0. Although there is only Ry + Rs left for R2 in either 
of the two limits, the latter (sharp interface) limit imposes vanishing interfacial normal 
velocity as the additional condition. This would certainly lead to a flow whose total rate 
of dissipation is larger than that obtained from minimizing the same functional without 
the additional constraint. 

We want to point out that while mathematically the sharp interface limit may be 
simply obtained by excluding Rd and Rr from R2 and applying = at the interface 
instead, physically the two-phase interfacial dissipation may not always be negligible, 
especially when the interfacial region of partial miscibility has finite and non-negligible 
width such that structural relaxation may occur. 

It should be noted that dimensional analysis indicates that [M] = [T] [Length]'^, i.e., 
there is a length scale lo which links these two parameters through the relation M = TIq. 
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Physically it is plausible to assume that lo is determined by a combination of microscopic 

factors, such as fluid-fluid interaction, fluid-solid interaction, molecular organization of 

the fluids, and molecular structure of the wall. Hence M and F may be physically related 

in any given system. 

6. Concluding remarks 

We should point out the inadequacies in our present formulation. First, the free energy 
used to delineate the two fluids is a minimal model. It neglects, for example, the density 
variation that can be fairly significant in the interfacial region. Second, our model in its 
present form is only applicable to simple liquids. Complex fluids would require nontrivial 
extensions. Third, we have neglected the van der Waals interaction which is very impor- 
tant in understanding precursor films. These and other inadequacies represent tasks still 
to be pursued. The main purpose of this paper is to outline the framework of a general 
theory which can resolve the MCL problem in its simplest form. 

It is important to emphasize that while the partial slip in single-phase flows is generally 
small and quantitatively indistinguishable from no-slip, yet its significance is qualitatively 
much greater. First, it is clear from the above that the GNBC goes hand-in-hand with 
the NBC in single-phase flows, and that it is incompatible with the no-slip boundary 
condition even in single-phase flows. The latter is clear from the power-law partial slip 
which extends mesoscopic distances into single-phase flow regimes. Second, even if the 
slip is small, the fact that slip exists means that its magnitude may be manipulated, i.e., it 
can be made larger or smaller. In particular, since the slip coefficient is a thermodynamic 
quantity, just as the viscosity, its magnitude should depend on molecular interactions 
and interface geometries (as well as the state variables such as temperature) (Barrat 
& Bocquet 1999b; Leger 2003; Granick, Zhu & Lee 2003; Zhu & Granick 2004; Neto 
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et al. 2005), a fact which can be used to advantage experimentally through nanoscale 
manipulations and environmental controls. For example, effective slip at nano-patterned 
surfaces has already been studied (Philip 1972; Lauga & Stone 2003; Cottin-Bizonne 
et al. 2003; Cottin-Bizonne et al. 2004; Priezjev, Darhuber & Troian 2005; Qian, Wang 
& Sheng 2005). In contrast, no-slip boundary condition is a clean-cut statement, with 
no room for adjustment or for physics considerations, only carries with it the burden 
of proof. Here the broad applicability of the no-slip boundary condition can not be 
considered as proof against slipping, as a very small amount of partial slip would clearly 
lead to similar results. It is therefore rather obvious that whereas slip/partial-slip can be 
derived from general principles and demonstrated through MD simulations, no-slip has 
yet to be proved with similar generalities. In fact, at present slip is already a subject with 
some fairly extensive literature. We refer to the reviews by Granick, Zhu & Lee (2003) 
and by Neto et al. (2005) for a more complete list of references. 

In closing, we note that in the case of the moving contact line, complete slip occurs 
in the linear regime (i.e., /3 is a constant). This is in contrast to the view that complete 
slip can only occur when "interface fracturing" occurs, i.e., in the fully nonlinear regime 
(sometimes also denoted as "super slipping" or "threshold slipping"), corresponding to a 
stress-dependent (3 (see e.g. Thompson & Troian 1997). It is rather likely that a statistical 
mechanical study of the slip coefficient can indeed produce a threshold behavior, e.g., /3 
approaches zero as the shear stress exceeds a certain threshold. However, in the linear 
regime this is not the case. The difference between the MCL complete slip and the 
interface fracturing is that in the case of MCL, the (complete) slip velocity can be very 
small and the relevant shear stress can be low. The fact that complete slip can occur in 
this case is due to the localized nature of the uncompensated Young stress, in addition 
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to the tangential viscous stress. Thus we can have complete slip in both the linear and 

nonlinear regimes, with different underlying physics. 
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Z0138. 



Appendix A. The principle of minimum energy dissipation 

To outline the principle of minimum energy dissipation (Onsager 1931a and 1931b), 
consider a system described by one single variable a, governed by the overdamped 
Langevin equation 

dF{a) 



7a : 



da 



+at), 



(Al) 



where 7 is the frictional coefficient, a is the rate of change of a, F{a) is the free energy, 
and is a white noise satisfying {({t)({t')) — 2jkBT5{t — t'), with ks denoting the 
Boltzmann constant and T the temperature. The probability of finding the system in 
the state described by a is a function of time, denoted by P{a, t) and governed by the 
Fokker-Planck equation 
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where D is the diffusion coefficient satisfying the Einstein relation 7D = fcsT. It is clear 
that the Boltzmann distribution Peq{a) oc exp [—F{a)/kBT] is a stationary solution of 
the Fokker-Planck equation (A 2). It can be shown (Langer 1968) that the transition 
probability from a at t to a' a,t t + At, i.e., P2{a', t + At; a, t), is given by 



P2{a',t + At;a,t) 
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for a' in the vicinity of a and short time interval At. The most probable transition occurs 



between a and a' is the one which minimizes 



2 



7.2 ^ . 



At, (A 4) 



Here a = (a' — a)/ At and the minimum of A is taken with respect to a' , or equivalently 
d, for prescribed a. The Euler-Lagrange equation for minimizing A is thus 

as expected from the Langevin equation (Al). Equation (A 5) is actually the simplest, 
one- variable version of the linear relation in equation (3.1) for rates and forces, and the 
function 'ya'^/2 + {dF/da)a in A is the corresponding one- variable version of the function 
$(d, Q.) + F{a, a) in Onsager's variational principle, stated by equations (3.2), (3.3), and 
(3.4). From the above discussion, it is clear that (1) the variation of A should be taken 
with respect to the rate a, for prescribed state variable a, (2) the minimum dissipation 
principle implies the balance of dissipative force and the force derived from free energy 
(see equation (A 5)), and (3) the minimum dissipation principle yields the most probable 
course of a dissipative process, provided the displacement from the equilibrium is small 
(Onsager 1931b, Onsager & Machlup 1953). 



Appendix B. The dissipation function in steady state 

In a steady state with d^/dt = 0, the total rate of energy dissipation R2 = 2$ in 
equations (3.21) and (3.22) can be expressed as a functional of v(r) only. This is because 
for prescribed phase field (j), the rate J in the bulk can be determined from cj) = \ ■ \J(j) 
through equations (3.28) and (3.29) while the rate (j) = Vrdr4> at the solid surface is 
already given in terms of v. By expressing R2 as a functional of v(r), we can obtain a 
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form of the dissipation function that is useful to the variational analysis of the sharp 

interface limit and the diffuse interface limit. 

Prom equation (3.29) with d(j)/dt — 0, we can formally express /x(r) by 



(Bl) 



where G(r, r') is the Green function for the Laplacian operator satisfying the boundary 
condition = at the solid surface. From in equation (3.19) and J = —MVfi 

(equation (3.28)), we have 



Ra = j dr [M(V/x)2] = j dr [-M/xV^] 



(B2) 



where the integration by parts has been used with = at the solid surface. Substi- 
tuting equation (B 1) into (B 2) yields 



[v] = j dTidT2 



1 



--G(ri,r2) (v.V0)r, (v • V.^), 



(B3) 



Meanwhile, substituting </> = Vrdr<p into Rr in equation (3.20), we obtain 

R^ [v] = j dS 



^ {Vrdr(t)f 



(B4) 



Combining equations (3.9), (3.10), (B3), and (B4) with R2 defined in equations (3.21) 
and (3.22), we obtain 

i?2[v] = R,[^]+R,[v] + Ra[v]+RrM 

+ 



J dridr2 



1 



-G(ri,r2)(v.V</.)r, (v • V^), 



JdS 



(B5) 



for the total rate of dissipation in two-phase flows, which is the sum of Ry due to viscosity, 
Rs due to slipping, R^ due to diffusion in the bulk, and Rr due to relaxation at the surface. 

In accordance with the principle of minimum energy dissipation (Oiisager 1931a and 
1931b), for steady state the equation(s) of motion can be derived from minimizing the 
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dissipation-function $ (= |i?2 here) with respect to the rates. Here we note that in 
steady state, the rates J in the bulk and ^ at the sohd surface are already determined 
by V for prescribed and i?2[v] in equation (B5) displays a symmetric, quadratic form 
as a function(al) of v, in accordance with the reciprocal relations, i.e., pij = pji for the 
coefficients pij in ^{a,a) = ^J2i,jPij^i'^j (equation (3.3)). The variation of the total 
dissipation R2 [v] should be taken with respect to v only, as it is the only rate. 

Based on equation (B 5), minimizing R2 with respect to v in the bulk yields the Stokes 
equation (3.31), while minimizing R2 with respect to tangential fluid velocity Vr at the 
solid surface yields the GNBC (3.32). That is, consider a variation of the velocity field 
v(r) — > v(r)-|-5v(r). The associated changes in Ry and Rs are already given by equations 
(3.11) and (3.12), and those in ii<i[v] and Rr[v] are given by 

and 

5Rr = j dS ^ {Vrdr(f>) 5Vrdr(t) = -2 j dS [L{(f>)dr(f>SVr] , (B 7) 

where equations (3.29) and (3.30) have been used. Combining equations (3.11), (3.12), 
(3.13), (B6), and (B7), we obtain the Euler-Lagrange equations 

-2r]dj {djVi + diVj) - diU - 2iidi(l) = (B 8) 

in the bulk and 

277 {dnVr + drVn) + 2(3vfP - 2L{4>)dr<^ = (B 9) 

at the surface. Note that equation (B8) is identical to the Stokes equation (3.31) with 
djVj = and a = —2p, and equation (B9) reduces to the GNBC (3.32). 
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